*************************************************************
*** Analyze the Monte Carlo experiments (SLX data-generating process) for the WWW JoP project
***
*** Created: 9-19-14
*** Modified: 9-5-19
***
*************************************************************

*** Open the data containing the SLX simulations
use "Experiments\SLX\SLX.dta", clear

*** Merge in the true SLX effects dataset
sort theta_1 theta_2
merge theta_1 theta_2 using "Experiments\SLX\SLX--True Effects.dta"
drop _merge

preserve
	duplicates drop theta_1 theta_2, force
	gen id = _n
	sort theta_1 theta_2
	tempfile id
	save `id', replace 
restore

sort theta_1 theta_2
merge theta_1 theta_2 using `id'
drop _merge

qui sum id
local end = r(max)

tempname slx
postfile `slx' str20 model modelid theta1 theta2 reject_p0 /*
*/	rec_teffect_x1 rec_deffect_x1 rec_ieffect_x1 rec_feffect_x1 rec_heffect_x1 rec_deffect_o0_x1 rec_teffect_o1_x1 rec_deffect_o1_x1 rec_ieffect_o1_x1 rec_teffect_o2_x1 rec_deffect_o2_x1 rec_ieffect_o2_x1 rec_teffect_o3_x1 rec_deffect_o3_x1 rec_ieffect_o3_x1 /*
*/	rec_teffect_x2 rec_deffect_x2 rec_ieffect_x2 rec_feffect_x2 rec_heffect_x2 rec_deffect_o0_x2 rec_teffect_o1_x2 rec_deffect_o1_x2 rec_ieffect_o1_x2 rec_teffect_o2_x2 rec_deffect_o2_x2 rec_ieffect_o2_x2 rec_teffect_o3_x2 rec_deffect_o3_x2 rec_ieffect_o3_x2 /*
*/	avg_teffect_x1 avg_deffect_x1 avg_ieffect_x1 avg_feffect_x1 avg_heffect_x1 avg_deffect_o0_x1 avg_teffect_o1_x1 avg_deffect_o1_x1 avg_ieffect_o1_x1 avg_teffect_o2_x1 avg_deffect_o2_x1 avg_ieffect_o2_x1 avg_teffect_o3_x1 avg_deffect_o3_x1 avg_ieffect_o3_x1 /*
*/	avg_teffect_x2 avg_deffect_x2 avg_ieffect_x2 avg_feffect_x2 avg_heffect_x2 avg_deffect_o0_x2 avg_teffect_o1_x2 avg_deffect_o1_x2 avg_ieffect_o1_x2 avg_teffect_o2_x2 avg_deffect_o2_x2 avg_ieffect_o2_x2 avg_teffect_o3_x2 avg_deffect_o3_x2 avg_ieffect_o3_x2 /*
*/	true_teffect_x1 true_deffect_x1 true_ieffect_x1 true_feffect_x1 true_heffect_x1 true_deffect_o0_x1 true_teffect_o1_x1 true_deffect_o1_x1 true_ieffect_o1_x1 true_teffect_o2_x1 true_deffect_o2_x1 true_ieffect_o2_x1 true_teffect_o3_x1 true_deffect_o3_x1 true_ieffect_o3_x1 /*
*/	true_teffect_x2 true_deffect_x2 true_ieffect_x2 true_feffect_x2 true_heffect_x2 true_deffect_o0_x2 true_teffect_o1_x2 true_deffect_o1_x2 true_ieffect_o1_x2 true_teffect_o2_x2 true_deffect_o2_x2 true_ieffect_o2_x2 true_teffect_o3_x2 true_deffect_o3_x2 true_ieffect_o3_x2 /*
*/	using "Experiments\Analyze\SLX\Analyze SLX Monte Carlo Experiments.dta", replace 

*** Model 1: underspecified SLX
preserve
	local m = 1
	
	qui foreach x of numlist 1 2 {
		tempvar sig_teffect_x`x' sig_deffect_x`x' sig_deffect_o0_x`x' 
		gen `sig_teffect_x`x'' = cond((teffect_x`x'_m`m' - 1.96*teffect_x`x'_se_m`m' < true_teffect_x`x') & (teffect_x`x'_m`m' + 1.96*teffect_x`x'_se_m`m' > true_teffect_x`x'), 1, 0)
		gen `sig_deffect_x`x'' = cond((deffect_x`x'_m`m' - 1.96*deffect_x`x'_se_m`m' < true_deffect_x`x') & (deffect_x`x'_m`m' + 1.96*deffect_x`x'_se_m`m' > true_deffect_x`x'), 1, 0)
		gen `sig_deffect_o0_x`x'' = cond((deffect_o0_x`x'_m`m' - 1.96*deffect_o0_x`x'_se_m`m' < true_deffect_o0_x`x') & (deffect_o0_x`x'_m`m' + 1.96*deffect_o0_x`x'_se_m`m' > true_deffect_o0_x`x'), 1, 0)		
	}

	qui foreach n of numlist 1(1)`end' {

		foreach v of varlist true_teffect_x1 - true_ieffect_o3_x2 {
			sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}	
	
		foreach x of numlist 1 2 {	
			sum theta_`x' if id == `n', meanonly
			local theta_`x' = r(mean)
			
			foreach v in teffect_x`x' deffect_x`x' deffect_o0_x`x' {
				sum `v'_m`m' if id == `n', meanonly
				local avg_`v' = r(mean)
			
				sum `sig_`v'' if id == `n', meanonly
				local rec_`v' = 100 * r(mean)
			}
		}		
		
		nois di _newline(2) "******************************************************"
		nois di "Theta_1 = `theta_1'; Theta_2 = `theta_2'"
		nois di "% of cases where Model 1 (underspecified) recovers the true total effect for x1 = " `rec_teffect_x1'
		nois di "% of cases where Model 1 (underspecified) recovers the true direct effect for x1 = " `rec_deffect_x1'
		nois di "% of cases where Model 1 (underspecified) recovers the true zero-order direct effect for x1 = " `rec_deffect_o0_x1'
		nois di _newline(1)
		nois di "% of cases where Model 1 (underspecified) recovers the true total effect for x2 = " `rec_teffect_x2'
		nois di "% of cases where Model 1 (underspecified) recovers the true direct effect for x2 = " `rec_deffect_x2'
		nois di "% of cases where Model 1 (underspecified) recovers the true zero-order direct effect for x2 = " `rec_deffect_o0_x2'
		nois di _newline(2)
		
		post `slx' ("SLX1") (`n') (`theta_1') (`theta_2') (.) /*		
		*/	(`rec_teffect_x1') (`rec_deffect_x1') (.) (.) (.) (`rec_deffect_o0_x1') (.) (.) (.) (.) (.) (.) (.) (.) (.) /*	
		*/	(`rec_teffect_x2') (`rec_deffect_x2') (.) (.) (.) (`rec_deffect_o0_x2') (.) (.) (.) (.) (.) (.) (.) (.) (.) /*	
		*/	(`avg_teffect_x1') (`avg_deffect_x1') (.) (.) (.) (`avg_deffect_o0_x1') (.) (.) (.) (.) (.) (.) (.) (.) (.) /*	
		*/	(`avg_teffect_x2') (`avg_deffect_x2') (.) (.) (.) (`avg_deffect_o0_x2') (.) (.) (.) (.) (.) (.) (.) (.) (.) /*	
		*/	(`true_teffect_x1') (`true_deffect_x1') (`true_ieffect_x1') (`true_feffect_x1') (`true_heffect_x1') (`true_deffect_o0_x1') (`true_teffect_o1_x1') (`true_deffect_o1_x1') (`true_ieffect_o1_x1') (`true_teffect_o2_x1') (`true_deffect_o2_x1') (`true_ieffect_o2_x1') (`true_teffect_o3_x1') (`true_deffect_o3_x1') (`true_ieffect_o3_x1') /*
		*/	(`true_teffect_x2') (`true_deffect_x2') (`true_ieffect_x2') (`true_feffect_x2') (`true_heffect_x2') (`true_deffect_o0_x2') (`true_teffect_o1_x2') (`true_deffect_o1_x2') (`true_ieffect_o1_x2') (`true_teffect_o2_x2') (`true_deffect_o2_x2') (`true_ieffect_o2_x2') (`true_teffect_o3_x2') (`true_deffect_o3_x2') (`true_ieffect_o3_x2') 	
	}
restore

*** Model 2: standard SLX (correct model)
preserve
	local m = 2

	qui foreach x of numlist 1 2 {
		tempvar sig_teffect_x`x' sig_deffect_x`x' sig_ieffect_x`x' sig_deffect_o0_x`x' sig_teffect_o1_x`x' sig_deffect_o1_x`x' sig_ieffect_o1_x`x'
		gen `sig_teffect_x`x'' = cond((teffect_x`x'_m`m' - 1.96*teffect_x`x'_se_m`m' < true_teffect_x`x') & (teffect_x`x'_m`m' + 1.96*teffect_x`x'_se_m`m' > true_teffect_x`x'), 1, 0)
		gen `sig_deffect_x`x'' = cond((deffect_x`x'_m`m' - 1.96*deffect_x`x'_se_m`m' < true_deffect_x`x') & (deffect_x`x'_m`m' + 1.96*deffect_x`x'_se_m`m' > true_deffect_x`x'), 1, 0)	
		gen `sig_ieffect_x`x'' = cond((ieffect_x`x'_m`m' - 1.96*ieffect_x`x'_se_m`m' < true_ieffect_x`x') & (ieffect_x`x'_m`m' + 1.96*ieffect_x`x'_se_m`m' > true_ieffect_x`x'), 1, 0)
		
		gen `sig_deffect_o0_x`x'' = cond((deffect_o0_x`x'_m`m' - 1.96*deffect_o0_x`x'_se_m`m' < true_deffect_o0_x`x') & (deffect_o0_x`x'_m`m' + 1.96*deffect_o0_x`x'_se_m`m' > true_deffect_o0_x`x'), 1, 0)		

		gen `sig_teffect_o1_x`x'' = cond((teffect_o1_x`x'_m`m' - 1.96*teffect_o1_x`x'_se_m`m' < true_teffect_o1_x`x') & (teffect_o1_x`x'_m`m' + 1.96*teffect_o1_x`x'_se_m`m' > true_teffect_o1_x`x'), 1, 0)
		gen `sig_deffect_o1_x`x'' = cond((deffect_o1_x`x'_m`m' - 1.96*deffect_o1_x`x'_se_m`m' < true_deffect_o1_x`x') & (deffect_o1_x`x'_m`m' + 1.96*deffect_o1_x`x'_se_m`m' > true_deffect_o1_x`x'), 1, 0)	
		gen `sig_ieffect_o1_x`x'' = cond((ieffect_o1_x`x'_m`m' - 1.96*ieffect_o1_x`x'_se_m`m' < true_ieffect_o1_x`x') & (ieffect_o1_x`x'_m`m' + 1.96*ieffect_o1_x`x'_se_m`m' > true_ieffect_o1_x`x'), 1, 0)		
	}

	qui foreach n of numlist 1(1)`end' {

		foreach v of varlist true_teffect_x1 - true_ieffect_o3_x2 {
			sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}	
	
		foreach x of numlist 1 2 {	
			sum theta_`x' if id == `n', meanonly
			local theta_`x' = r(mean)
			
			foreach v in teffect_x`x' deffect_x`x' ieffect_x`x' deffect_o0_x`x' teffect_o1_x`x' deffect_o1_x`x' ieffect_o1_x`x' {
				sum `v'_m`m' if id == `n', meanonly
				local avg_`v' = r(mean)
			
				sum `sig_`v'' if id == `n', meanonly
				local rec_`v' = 100 * r(mean)
			}
		}		
		
		nois di _newline(2) "******************************************************"
		nois di "Theta_1 = `theta_1'; Theta_2 = `theta_2'"
		nois di "% of cases where Model 2 (correct) recovers the true total effect for x1 = " `rec_teffect_x1'
		nois di "% of cases where Model 2 (correct) recovers the true indirect effect for x1 = " `rec_ieffect_x1'		
		nois di "% of cases where Model 2 (correct) recovers the true direct effect for x1 = " `rec_deffect_x1'
		nois di "% of cases where Model 2 (correct) recovers the true zero-order direct effect for x1 = " `rec_deffect_o0_x1'
		nois di "% of cases where Model 2 (correct) recovers the true first-order total effect for x1 = " `rec_teffect_o1_x1'
		nois di "% of cases where Model 2 (correct) recovers the true first-order indirect effect for x1 = " `rec_ieffect_o1_x1'		
		nois di _newline(1)
		nois di "% of cases where Model 2 (correct) recovers the true total effect for x2 = " `rec_teffect_x2'
		nois di "% of cases where Model 2 (correct) recovers the true indirect effect for x2 = " `rec_ieffect_x2'		
		nois di "% of cases where Model 2 (correct) recovers the true direct effect for x2 = " `rec_deffect_x2'
		nois di "% of cases where Model 2 (correct) recovers the true zero-order direct effect for x2 = " `rec_deffect_o0_x2'
		nois di "% of cases where Model 2 (correct) recovers the true first-order total effect for x2 = " `rec_teffect_o1_x2'
		nois di "% of cases where Model 2 (correct) recovers the true first-order indirect effect for x2 = " `rec_ieffect_o1_x2'		
		nois di _newline(2)
		
		post `slx' ("SLX2") (`n') (`theta_1') (`theta_2') (.) /*		
		*/	(`rec_teffect_x1') (`rec_deffect_x1') (`rec_ieffect_x1') (.) (.) (`rec_deffect_o0_x1') (`rec_teffect_o1_x1') (`rec_deffect_o1_x1') (`rec_ieffect_o1_x1') (.) (.) (.) (.) (.) (.) /*	
		*/	(`rec_teffect_x2') (`rec_deffect_x2') (`rec_ieffect_x2') (.) (.) (`rec_deffect_o0_x2') (`rec_teffect_o1_x2') (`rec_deffect_o1_x2') (`rec_ieffect_o1_x2') (.) (.) (.) (.) (.) (.) /*	
		*/	(`avg_teffect_x1') (`avg_deffect_x1') (`avg_ieffect_x1') (.) (.) (`avg_deffect_o0_x1') (`avg_teffect_o1_x1') (`avg_deffect_o1_x1') (`avg_ieffect_o1_x1') (.) (.) (.) (.) (.) (.) /*	
		*/	(`avg_teffect_x2') (`avg_deffect_x2') (`avg_ieffect_x2') (.) (.) (`avg_deffect_o0_x2') (`avg_teffect_o1_x2') (`avg_deffect_o1_x2') (`avg_ieffect_o1_x2') (.) (.) (.) (.) (.) (.) /*	
		*/	(`true_teffect_x1') (`true_deffect_x1') (`true_ieffect_x1') (`true_feffect_x1') (`true_heffect_x1') (`true_deffect_o0_x1') (`true_teffect_o1_x1') (`true_deffect_o1_x1') (`true_ieffect_o1_x1') (`true_teffect_o2_x1') (`true_deffect_o2_x1') (`true_ieffect_o2_x1') (`true_teffect_o3_x1') (`true_deffect_o3_x1') (`true_ieffect_o3_x1') /*
		*/	(`true_teffect_x2') (`true_deffect_x2') (`true_ieffect_x2') (`true_feffect_x2') (`true_heffect_x2') (`true_deffect_o0_x2') (`true_teffect_o1_x2') (`true_deffect_o1_x2') (`true_ieffect_o1_x2') (`true_teffect_o2_x2') (`true_deffect_o2_x2') (`true_ieffect_o2_x2') (`true_teffect_o3_x2') (`true_deffect_o3_x2') (`true_ieffect_o3_x2') 	
	}
restore



*** Model 3: standard SLX (overspecified model)
preserve
	local m = 3
		
	qui foreach x of numlist 1 2 {
		tempvar sig_teffect_x`x' sig_deffect_x`x' sig_ieffect_x`x' sig_feffect_x`x' sig_heffect_x`x' sig_deffect_o0_x`x' sig_teffect_o1_x`x' sig_deffect_o1_x`x' sig_ieffect_o1_x`x' sig_teffect_o2_x`x' sig_deffect_o2_x`x' sig_ieffect_o2_x`x'
		gen `sig_teffect_x`x'' = cond((teffect_x`x'_m`m' - 1.96*teffect_x`x'_se_m`m' < true_teffect_x`x') & (teffect_x`x'_m`m' + 1.96*teffect_x`x'_se_m`m' > true_teffect_x`x'), 1, 0)
		gen `sig_deffect_x`x'' = cond((deffect_x`x'_m`m' - 1.96*deffect_x`x'_se_m`m' < true_deffect_x`x') & (deffect_x`x'_m`m' + 1.96*deffect_x`x'_se_m`m' > true_deffect_x`x'), 1, 0)	
		gen `sig_ieffect_x`x'' = cond((ieffect_x`x'_m`m' - 1.96*ieffect_x`x'_se_m`m' < true_ieffect_x`x') & (ieffect_x`x'_m`m' + 1.96*ieffect_x`x'_se_m`m' > true_ieffect_x`x'), 1, 0)

		gen `sig_feffect_x`x'' = cond((feffect_x`x'_m`m' - 1.96*feffect_x`x'_se_m`m' < true_feffect_x`x') & (feffect_x`x'_m`m' + 1.96*feffect_x`x'_se_m`m' > true_feffect_x`x'), 1, 0)
		gen `sig_heffect_x`x'' = cond((heffect_x`x'_m`m' - 1.96*heffect_x`x'_se_m`m' < true_heffect_x`x') & (heffect_x`x'_m`m' + 1.96*heffect_x`x'_se_m`m' > true_heffect_x`x'), 1, 0)
		
		gen `sig_deffect_o0_x`x'' = cond((deffect_o0_x`x'_m`m' - 1.96*deffect_o0_x`x'_se_m`m' < true_deffect_o0_x`x') & (deffect_o0_x`x'_m`m' + 1.96*deffect_o0_x`x'_se_m`m' > true_deffect_o0_x`x'), 1, 0)		

		gen `sig_teffect_o1_x`x'' = cond((teffect_o1_x`x'_m`m' - 1.96*teffect_o1_x`x'_se_m`m' < true_teffect_o1_x`x') & (teffect_o1_x`x'_m`m' + 1.96*teffect_o1_x`x'_se_m`m' > true_teffect_o1_x`x'), 1, 0)
		gen `sig_deffect_o1_x`x'' = cond((deffect_o1_x`x'_m`m' - 1.96*deffect_o1_x`x'_se_m`m' < true_deffect_o1_x`x') & (deffect_o1_x`x'_m`m' + 1.96*deffect_o1_x`x'_se_m`m' > true_deffect_o1_x`x'), 1, 0)	
		gen `sig_ieffect_o1_x`x'' = cond((ieffect_o1_x`x'_m`m' - 1.96*ieffect_o1_x`x'_se_m`m' < true_ieffect_o1_x`x') & (ieffect_o1_x`x'_m`m' + 1.96*ieffect_o1_x`x'_se_m`m' > true_ieffect_o1_x`x'), 1, 0)		
		
		gen `sig_teffect_o2_x`x'' = cond((teffect_o2_x`x'_m`m' - 1.96*teffect_o2_x`x'_se_m`m' < true_teffect_o2_x`x') & (teffect_o2_x`x'_m`m' + 1.96*teffect_o2_x`x'_se_m`m' > true_teffect_o2_x`x'), 1, 0)
		gen `sig_deffect_o2_x`x'' = cond((deffect_o2_x`x'_m`m' - 1.96*deffect_o2_x`x'_se_m`m' < true_deffect_o2_x`x') & (deffect_o2_x`x'_m`m' + 1.96*deffect_o2_x`x'_se_m`m' > true_deffect_o2_x`x'), 1, 0)	
		gen `sig_ieffect_o2_x`x'' = cond((ieffect_o2_x`x'_m`m' - 1.96*ieffect_o2_x`x'_se_m`m' < true_ieffect_o2_x`x') & (ieffect_o2_x`x'_m`m' + 1.96*ieffect_o2_x`x'_se_m`m' > true_ieffect_o2_x`x'), 1, 0)				
	}

	qui foreach n of numlist 1(1)`end' {

		foreach v of varlist true_teffect_x1 - true_ieffect_o3_x2 {
			sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}	
	
		foreach x of numlist 1 2 {	
			sum theta_`x' if id == `n', meanonly
			local theta_`x' = r(mean)
			
			foreach v in teffect_x`x' deffect_x`x' ieffect_x`x' feffect_x`x' heffect_x`x' deffect_o0_x`x' teffect_o1_x`x' deffect_o1_x`x' ieffect_o1_x`x' teffect_o2_x`x' deffect_o2_x`x' ieffect_o2_x`x' {
				sum `v'_m`m' if id == `n', meanonly
				local avg_`v' = r(mean)
			
				sum `sig_`v'' if id == `n', meanonly
				local rec_`v' = 100 * r(mean)
			}
		}		
		
		nois di _newline(2) "******************************************************"
		nois di "Theta_1 = `theta_1'; Theta_2 = `theta_2'"
		nois di "% of cases where Model 3 (overspecified) recovers the true total effect for x1 = " `rec_teffect_x1'
		nois di "% of cases where Model 3 (overspecified) recovers the true indirect effect for x1 = " `rec_ieffect_x1'		
		nois di "% of cases where Model 3 (overspecified) recovers the true direct effect for x1 = " `rec_deffect_x1'
		nois di "% of cases where Model 3 (overspecified) recovers the true zero-order direct effect for x1 = " `rec_deffect_o0_x1'
		nois di "% of cases where Model 3 (overspecified) recovers the true first-order total effect for x1 = " `rec_teffect_o1_x1'
		nois di "% of cases where Model 3 (overspecified) recovers the true first-order indirect effect for x1 = " `rec_ieffect_o1_x1'		
		nois di "% of cases where Model 3 (overspecified) recovers the true second-order total effect for x1 = " `rec_teffect_o2_x1'
		nois di "% of cases where Model 3 (overspecified) recovers the true second-order direct effect for x1 = " `rec_deffect_o2_x1'		
		nois di "% of cases where Model 3 (overspecified) recovers the true second-order indirect effect for x1 = " `rec_ieffect_o2_x1'		
		
		nois di _newline(1)
		nois di "% of cases where Model 3 (overspecified) recovers the true total effect for x2 = " `rec_teffect_x2'
		nois di "% of cases where Model 3 (overspecified) recovers the true indirect effect for x2 = " `rec_ieffect_x2'		
		nois di "% of cases where Model 3 (overspecified) recovers the true direct effect for x2 = " `rec_deffect_x2'
		nois di "% of cases where Model 3 (overspecified) recovers the true zero-order direct effect for x2 = " `rec_deffect_o0_x2'
		nois di "% of cases where Model 3 (overspecified) recovers the true first-order total effect for x2 = " `rec_teffect_o1_x2'
		nois di "% of cases where Model 3 (overspecified) recovers the true first-order indirect effect for x2 = " `rec_ieffect_o1_x2'		
		nois di "% of cases where Model 3 (overspecified) recovers the true second-order total effect for x2 = " `rec_teffect_o2_x2'
		nois di "% of cases where Model 3 (overspecified) recovers the true second-order direct effect for x2 = " `rec_deffect_o2_x2'		
		nois di "% of cases where Model 3 (overspecified) recovers the true second-order indirect effect for x2 = " `rec_ieffect_o2_x2'				
		nois di _newline(2)
		
		post `slx' ("SLX3") (`n') (`theta_1') (`theta_2') (.) /*		
		*/	(`rec_teffect_x1') (`rec_deffect_x1') (`rec_ieffect_x1') (`rec_feffect_x1') (`rec_heffect_x1') (`rec_deffect_o0_x1') (`rec_teffect_o1_x1') (`rec_deffect_o1_x1') (`rec_ieffect_o1_x1') (`rec_teffect_o2_x1') (`rec_deffect_o2_x1') (`rec_ieffect_o2_x1') (.) (.) (.) /*	
		*/	(`rec_teffect_x2') (`rec_deffect_x2') (`rec_ieffect_x2') (`rec_feffect_x2') (`rec_heffect_x2') (`rec_deffect_o0_x2') (`rec_teffect_o1_x2') (`rec_deffect_o1_x2') (`rec_ieffect_o1_x2') (`rec_teffect_o2_x2') (`rec_deffect_o2_x2') (`rec_ieffect_o2_x2') (.) (.) (.) /*	
		*/	(`avg_teffect_x1') (`avg_deffect_x1') (`avg_ieffect_x1') (`avg_feffect_x1') (`avg_heffect_x1') (`avg_deffect_o0_x1') (`avg_teffect_o1_x1') (`avg_deffect_o1_x1') (`avg_ieffect_o1_x1') (`avg_teffect_o2_x1') (`avg_deffect_o2_x1') (`avg_ieffect_o2_x1') (.) (.) (.) /*	
		*/	(`avg_teffect_x2') (`avg_deffect_x2') (`avg_ieffect_x2') (`avg_feffect_x2') (`avg_heffect_x2') (`avg_deffect_o0_x2') (`avg_teffect_o1_x2') (`avg_deffect_o1_x2') (`avg_ieffect_o1_x2') (`avg_teffect_o2_x2') (`avg_deffect_o2_x2') (`avg_ieffect_o2_x2') (.) (.) (.) /*	
		*/	(`true_teffect_x1') (`true_deffect_x1') (`true_ieffect_x1') (`true_feffect_x1') (`true_heffect_x1') (`true_deffect_o0_x1') (`true_teffect_o1_x1') (`true_deffect_o1_x1') (`true_ieffect_o1_x1') (`true_teffect_o2_x1') (`true_deffect_o2_x1') (`true_ieffect_o2_x1') (`true_teffect_o3_x1') (`true_deffect_o3_x1') (`true_ieffect_o3_x1') /*
		*/	(`true_teffect_x2') (`true_deffect_x2') (`true_ieffect_x2') (`true_feffect_x2') (`true_heffect_x2') (`true_deffect_o0_x2') (`true_teffect_o1_x2') (`true_deffect_o1_x2') (`true_ieffect_o1_x2') (`true_teffect_o2_x2') (`true_deffect_o2_x2') (`true_ieffect_o2_x2') (`true_teffect_o3_x2') (`true_deffect_o3_x2') (`true_ieffect_o3_x2') 	
	}
restore



*** Model 4: SAR
preserve
	local m = 4

	qui foreach x of numlist 1 2 {
		tempvar sig_teffect_x`x' sig_deffect_x`x' sig_ieffect_x`x' sig_feffect_x`x' sig_heffect_x`x' sig_deffect_o0_x`x' sig_teffect_o1_x`x' sig_deffect_o1_x`x' sig_ieffect_o1_x`x' sig_teffect_o2_x`x' sig_deffect_o2_x`x' sig_ieffect_o2_x`x' sig_teffect_o3_x`x' sig_deffect_o3_x`x' sig_ieffect_o3_x`x'
		gen `sig_teffect_x`x'' = cond((teffect_x`x'_m`m' - 1.96*teffect_x`x'_se_m`m' < true_teffect_x`x') & (teffect_x`x'_m`m' + 1.96*teffect_x`x'_se_m`m' > true_teffect_x`x'), 1, 0)
		gen `sig_deffect_x`x'' = cond((deffect_x`x'_m`m' - 1.96*deffect_x`x'_se_m`m' < true_deffect_x`x') & (deffect_x`x'_m`m' + 1.96*deffect_x`x'_se_m`m' > true_deffect_x`x'), 1, 0)	
		gen `sig_ieffect_x`x'' = cond((ieffect_x`x'_m`m' - 1.96*ieffect_x`x'_se_m`m' < true_ieffect_x`x') & (ieffect_x`x'_m`m' + 1.96*ieffect_x`x'_se_m`m' > true_ieffect_x`x'), 1, 0)

		gen `sig_feffect_x`x'' = cond((feffect_x`x'_m`m' - 1.96*feffect_x`x'_se_m`m' < true_feffect_x`x') & (feffect_x`x'_m`m' + 1.96*feffect_x`x'_se_m`m' > true_feffect_x`x'), 1, 0)
		gen `sig_heffect_x`x'' = cond((heffect_x`x'_m`m' - 1.96*heffect_x`x'_se_m`m' < true_heffect_x`x') & (heffect_x`x'_m`m' + 1.96*heffect_x`x'_se_m`m' > true_heffect_x`x'), 1, 0)
		
		gen `sig_deffect_o0_x`x'' = cond((deffect_o0_x`x'_m`m' - 1.96*deffect_o0_x`x'_se_m`m' < true_deffect_o0_x`x') & (deffect_o0_x`x'_m`m' + 1.96*deffect_o0_x`x'_se_m`m' > true_deffect_o0_x`x'), 1, 0)		

		gen `sig_teffect_o1_x`x'' = cond((teffect_o1_x`x'_m`m' - 1.96*teffect_o1_x`x'_se_m`m' < true_teffect_o1_x`x') & (teffect_o1_x`x'_m`m' + 1.96*teffect_o1_x`x'_se_m`m' > true_teffect_o1_x`x'), 1, 0)
		gen `sig_deffect_o1_x`x'' = cond((deffect_o1_x`x'_m`m' - 1.96*deffect_o1_x`x'_se_m`m' < true_deffect_o1_x`x') & (deffect_o1_x`x'_m`m' + 1.96*deffect_o1_x`x'_se_m`m' > true_deffect_o1_x`x'), 1, 0)	
		gen `sig_ieffect_o1_x`x'' = cond((ieffect_o1_x`x'_m`m' - 1.96*ieffect_o1_x`x'_se_m`m' < true_ieffect_o1_x`x') & (ieffect_o1_x`x'_m`m' + 1.96*ieffect_o1_x`x'_se_m`m' > true_ieffect_o1_x`x'), 1, 0)		
		
		gen `sig_teffect_o2_x`x'' = cond((teffect_o2_x`x'_m`m' - 1.96*teffect_o2_x`x'_se_m`m' < true_teffect_o2_x`x') & (teffect_o2_x`x'_m`m' + 1.96*teffect_o2_x`x'_se_m`m' > true_teffect_o2_x`x'), 1, 0)
		gen `sig_deffect_o2_x`x'' = cond((deffect_o2_x`x'_m`m' - 1.96*deffect_o2_x`x'_se_m`m' < true_deffect_o2_x`x') & (deffect_o2_x`x'_m`m' + 1.96*deffect_o2_x`x'_se_m`m' > true_deffect_o2_x`x'), 1, 0)	
		gen `sig_ieffect_o2_x`x'' = cond((ieffect_o2_x`x'_m`m' - 1.96*ieffect_o2_x`x'_se_m`m' < true_ieffect_o2_x`x') & (ieffect_o2_x`x'_m`m' + 1.96*ieffect_o2_x`x'_se_m`m' > true_ieffect_o2_x`x'), 1, 0)				
		
		gen `sig_teffect_o3_x`x'' = cond((teffect_o3_x`x'_m`m' - 1.96*teffect_o3_x`x'_se_m`m' < true_teffect_o3_x`x') & (teffect_o3_x`x'_m`m' + 1.96*teffect_o3_x`x'_se_m`m' > true_teffect_o3_x`x'), 1, 0)
		gen `sig_deffect_o3_x`x'' = cond((deffect_o3_x`x'_m`m' - 1.96*deffect_o3_x`x'_se_m`m' < true_deffect_o3_x`x') & (deffect_o3_x`x'_m`m' + 1.96*deffect_o3_x`x'_se_m`m' > true_deffect_o3_x`x'), 1, 0)	
		gen `sig_ieffect_o3_x`x'' = cond((ieffect_o3_x`x'_m`m' - 1.96*ieffect_o3_x`x'_se_m`m' < true_ieffect_o3_x`x') & (ieffect_o3_x`x'_m`m' + 1.96*ieffect_o3_x`x'_se_m`m' > true_ieffect_o3_x`x'), 1, 0)						
	}

	tempvar sig_rho
	gen `sig_rho' = cond((rho_m4 - 1.96*rho_se_m4 < 0) & (rho_m4 + 1.96*rho_se_m4 > 0), 0, 1)
	
	qui foreach n of numlist 1(1)`end' {

		foreach v of varlist true_teffect_x1 - true_ieffect_o3_x2 {
			sum `v' if id == `n', meanonly
			local `v' = r(mean)
		}	
	
		foreach x of numlist 1 2 {	
			sum theta_`x' if id == `n', meanonly
			local theta_`x' = r(mean)
			
			foreach v in teffect_x`x' deffect_x`x' ieffect_x`x' feffect_x`x' heffect_x`x' deffect_o0_x`x' teffect_o1_x`x' deffect_o1_x`x' ieffect_o1_x`x' teffect_o2_x`x' deffect_o2_x`x' ieffect_o2_x`x' teffect_o3_x`x' deffect_o3_x`x' ieffect_o3_x`x' {
				sum `v'_m`m' if id == `n', meanonly
				local avg_`v' = r(mean)
			
				sum `sig_`v'' if id == `n', meanonly
				local rec_`v' = 100 * r(mean)
			}
		}		
		
		sum `sig_rho' if id == `n', meanonly
		local reject_p0 = 100 * r(mean)
		
		nois di _newline(2) "******************************************************"
		nois di "Theta_1 = `theta_1'; Theta_2 = `theta_2'"
		nois di "% of cases where Model 4 (SAR) recovers the true total effect for x1 = " `rec_teffect_x1'
		nois di "% of cases where Model 4 (SAR) recovers the true indirect effect for x1 = " `rec_ieffect_x1'		
		nois di "% of cases where Model 4 (SAR) recovers the true direct effect for x1 = " `rec_deffect_x1'
		nois di "% of cases where Model 4 (SAR) recovers the true zero-order direct effect for x1 = " `rec_deffect_o0_x1'
		nois di "% of cases where Model 4 (SAR) recovers the true first-order total effect for x1 = " `rec_teffect_o1_x1'
		nois di "% of cases where Model 4 (SAR) recovers the true first-order indirect effect for x1 = " `rec_ieffect_o1_x1'		
		nois di "% of cases where Model 4 (SAR) recovers the true second-order total effect for x1 = " `rec_teffect_o2_x1'
		nois di "% of cases where Model 4 (SAR) recovers the true second-order direct effect for x1 = " `rec_deffect_o2_x1'		
		nois di "% of cases where Model 4 (SAR) recovers the true second-order indirect effect for x1 = " `rec_ieffect_o2_x1'		
		nois di "% of cases where Model 4 (SAR) recovers the true feedback effects for x1 = " `rec_feffect_x1'		
		nois di "% of cases where Model 4 (SAR) recovers the true higher-order effects for x1 = " `rec_heffect_x1'				
		
		nois di _newline(1)
		nois di "% of cases where Model 4 (SAR) recovers the true total effect for x2 = " `rec_teffect_x2'
		nois di "% of cases where Model 4 (SAR) recovers the true indirect effect for x2 = " `rec_ieffect_x2'		
		nois di "% of cases where Model 4 (SAR) recovers the true direct effect for x2 = " `rec_deffect_x2'
		nois di "% of cases where Model 4 (SAR) recovers the true zero-order direct effect for x2 = " `rec_deffect_o0_x2'
		nois di "% of cases where Model 4 (SAR) recovers the true first-order total effect for x2 = " `rec_teffect_o1_x2'
		nois di "% of cases where Model 4 (SAR) recovers the true first-order indirect effect for x2 = " `rec_ieffect_o1_x2'		
		nois di "% of cases where Model 4 (SAR) recovers the true second-order total effect for x2 = " `rec_teffect_o2_x2'
		nois di "% of cases where Model 4 (SAR) recovers the true second-order direct effect for x2 = " `rec_deffect_o2_x2'		
		nois di "% of cases where Model 4 (SAR) recovers the true second-order indirect effect for x2 = " `rec_ieffect_o2_x2'				
		nois di "% of cases where Model 4 (SAR) recovers the true feedback effects for x2 = " `rec_feffect_x2'		
		nois di "% of cases where Model 4 (SAR) recovers the true higher-order effects for x2 = " `rec_heffect_x2'						
		nois di _newline(2)
		
		post `slx' ("SAR") (`n') (`theta_1') (`theta_2') (`reject_p0') /*		
		*/	(`rec_teffect_x1') (`rec_deffect_x1') (`rec_ieffect_x1') (`rec_feffect_x1') (`rec_heffect_x1') (`rec_deffect_o0_x1') (`rec_teffect_o1_x1') (`rec_deffect_o1_x1') (`rec_ieffect_o1_x1') (`rec_teffect_o2_x1') (`rec_deffect_o2_x1') (`rec_ieffect_o2_x1') (`rec_teffect_o3_x1') (`rec_deffect_o3_x1') (`rec_ieffect_o3_x1') /*	
		*/	(`rec_teffect_x2') (`rec_deffect_x2') (`rec_ieffect_x2') (`rec_feffect_x2') (`rec_heffect_x2') (`rec_deffect_o0_x2') (`rec_teffect_o1_x2') (`rec_deffect_o1_x2') (`rec_ieffect_o1_x2') (`rec_teffect_o2_x2') (`rec_deffect_o2_x2') (`rec_ieffect_o2_x2') (`rec_teffect_o3_x2') (`rec_deffect_o3_x2') (`rec_ieffect_o3_x2') /*		
		*/	(`avg_teffect_x1') (`avg_deffect_x1') (`avg_ieffect_x1') (`avg_feffect_x1') (`avg_heffect_x1') (`avg_deffect_o0_x1') (`avg_teffect_o1_x1') (`avg_deffect_o1_x1') (`avg_ieffect_o1_x1') (`avg_teffect_o2_x1') (`avg_deffect_o2_x1') (`avg_ieffect_o2_x1') (`avg_teffect_o3_x1') (`avg_deffect_o3_x1') (`avg_ieffect_o3_x1') /*	
		*/	(`avg_teffect_x2') (`avg_deffect_x2') (`avg_ieffect_x2') (`avg_feffect_x2') (`avg_heffect_x2') (`avg_deffect_o0_x2') (`avg_teffect_o1_x2') (`avg_deffect_o1_x2') (`avg_ieffect_o1_x2') (`avg_teffect_o2_x2') (`avg_deffect_o2_x2') (`avg_ieffect_o2_x2') (`avg_teffect_o3_x2') (`avg_deffect_o3_x2') (`avg_ieffect_o3_x2') /*	
		*/	(`true_teffect_x1') (`true_deffect_x1') (`true_ieffect_x1') (`true_feffect_x1') (`true_heffect_x1') (`true_deffect_o0_x1') (`true_teffect_o1_x1') (`true_deffect_o1_x1') (`true_ieffect_o1_x1') (`true_teffect_o2_x1') (`true_deffect_o2_x1') (`true_ieffect_o2_x1') (`true_teffect_o3_x1') (`true_deffect_o3_x1') (`true_ieffect_o3_x1') /*
		*/	(`true_teffect_x2') (`true_deffect_x2') (`true_ieffect_x2') (`true_feffect_x2') (`true_heffect_x2') (`true_deffect_o0_x2') (`true_teffect_o1_x2') (`true_deffect_o1_x2') (`true_ieffect_o1_x2') (`true_teffect_o2_x2') (`true_deffect_o2_x2') (`true_ieffect_o2_x2') (`true_teffect_o3_x2') (`true_deffect_o3_x2') (`true_ieffect_o3_x2') 	
	}
restore



postclose `slx'

use "Experiments\Analyze\SLX\Analyze SLX Monte Carlo Experiments.dta", clear


***************************** Direct Effects *****************************
list theta1 theta2 rec_deffect_o0_x1 rec_deffect_o1_x1 rec_deffect_o2_x1 rec_deffect_o0_x2 rec_deffect_o1_x2 rec_deffect_o2_x2 if model == "SAR"

***************************** Indirect Effects *****************************
list theta1 theta2 rec_ieffect_o1_x1 rec_ieffect_o2_x1 rec_ieffect_o1_x2 rec_ieffect_o2_x2 if model == "SAR"

***************************** Higher-Order Effects *****************************
list theta1 theta2 rec_feffect_x1 rec_heffect_x1 rec_feffect_x2 rec_heffect_x2 if model == "SAR"
